import itertools

# fct it's just for computation of \gamma(sigma_i)(I) for I=(i,J,i+1,K)
def fct(J,I):
    return([[I[0]+1,*I[1:I[0]],J[0],*I[I[0]:]],[*I[:I[0]+1],J[0],*I[I[0]+1:]]])

# calcplus(i,I) takes as input a positive integer "i" and a sequence "I=[i_1, ...,i_l]". Then compute \gamma(sigma_i)(i_1, ..., i_l)=coef_1[I1]+coef_2[I_2]+...+coef_m[I_m] and return the result in the form [[coef_1,I_1],[coef_2,I2], ..., [coef_m,I_m]]
def calcplus (i,I):
    k=-1
    l=-1
    j=0
    while (k+1)*(l+1)==0 and j<len(I):
        if I[j]==i :
            k=j
        if I[j]==i+1:
            l=j
        j+=1    
    if l==-1:
        if k==-1:
            return([[1]+I])
        return([[1,*I[:k],i+1,*I[k+1:]]])            
    if k==-1:
        I0=[*I[:l],i,*I[l+1:]]
        I1=[*I[:l],i,*I[l:]]
        if l>0:
            I2=[*I[:l+1],i,*I[l+1:]]
            return([[1]+I0,[1]+I1,[-1]+I2])
        return([[1]+I0,[1]+I1])
    if k>0:
        if k<l: 
            return([[1,*I[:k],i+1,*I[k+1:l],i,*I[l+1:]]])
        return([[1,*I[:l],i,*I[l+1:k],i+1,*I[k+1:]]])
    J=I[1:l]
    J.reverse()
    L=[[2,i,*I[l:]]]
    while J!=[]:
        L=[fct(J,K)[j] for K in L for j in range(0,2)]
        J=J[1:]
    for K in L:
        K[0]=(-1)**(K[0]+1)
    return(L) 

# calcminus(i,I) takes as input an integer "i" and a sequence "I=[i_1, ...,i_l]". Then it compute \gamma(sigma_i^{-1})(i_1, ..., i_l)=coef_1[I1]+coef_2[I_2]+...+coef_m[I_m] and return the result in the form [[coef_1,I_1],[coef_2,I2], ..., [coef_m,I_m]]
def calcminus (i,I):
    k=-1
    l=-1
    j=0
    while (k+1)*(l+1)==0 and j<len(I):
        if I[j]==i :
            k=j
        if I[j]==i+1:
            l=j
        j+=1    
    if k==-1:
        if l==-1:
            return([[1]+I])
        return([[1,*I[:l],i,*I[l+1:]]])            
    if l==-1:
        I0=[*I[:k],i+1,*I[k+1:]]
        I1=[*I[:k+1],i+1,*I[k+1:]]
        if k>0:
            I2=[*I[:k],i+1,*I[k:]]
            return([[1]+I0,[1]+I1,[-1]+I2])
        return([[1]+I0,[1]+I1])
    if k>0:
        if k<l: 
            return([[1,*I[:k],i+1,*I[k+1:l],i,*I[l+1:]]])
        return([[1,*I[:l],i,*I[l+1:k],i+1,*I[k+1:]]])
    J=I[1:l]
    J.reverse()
    L=[[2,i,*I[l:]]]
    while J!=[]:
        L=[fct(J,K)[j] for K in L for j in range(0,2)]
        J=J[1:]
    for K in L:
        K[0]=(-1)**(K[0]+1)
    return(L) 

# proj_calcgene (p,i,C) takes as input a positive integer "p", a relative integer "i" and a sequence of sequence "C=[[coef_1,I_1], ..., [coef_m,I_m]]". Then it computate gamma_p(sigma_{|i|}^{sign(i)})([coef_1,I_1], ..., [coef_m,I_m])=coef'_1[J_1]+...+coef'_l[J_l] and gives the result in the form [[coef'_1,J_1],[coef'_2,J_2],...,[coef'_l,J_l]]. (recall that gamma_p is the projection of the rpz on the subspace generated by commutator of weight lower or equal than p.)
def proj_calcgene (p,i,C):
    if i>0:
        return([[Y[0]*I[0]]+Y[1:] for I in C for Y in calcplus(i,I[1:]) if len(I)<=p+1 if len(Y)<=p+1])
    if i<0:
        return([[Y[0]*I[0]]+Y[1:] for I in C for Y in calcminus(-i,I[1:]) if len(I)<=p+1 if len(Y)<=p+1])
    return([[]])

# proj_calcgene (p,i,C) take as input a positive integer "p", a sequence of relatives integer "L=[l_1, ..., l_k" and a sequence of sequence "C=[[coef_1,I_1], ..., [coef_m,I_m]]". Then it computate gamma_p(sigma_{|l_1|}^{sign(l_1)}sigma_{|l_2|}^{sign(l_2)}...sigma_{|l_k|}^{sign(l_k)})([coef_1,I_1], ..., [coef_m,I_m])=coef'_1[J_1]+...+coef'_l[J_l] and gives the result in the form [[coef'_1,J_1],[coef'_2,J_2],...,[coef'_l,J_l]]. (recall that gamma_p is the projection of the rpz on the subspace generated by commutator of weight lower or equal than p.)

def proj_calcrpz(p,L,C):
    L.reverse()
    l=len(L)
    for i in range(0,l):
        C=proj_calcgene(p,L[i],C)
        j=0
        while j<len(C):
            k=j+1
            while k<len(C):
                if C[j][1:]==C[k][1:]:
                    C[j][0]=C[j][0]+C[k][0]
                    C.pop(k)
                    k-=1
                k+=1
            if C[j][0]==0:
                C.pop(j)
                j-=1
            j+=1
    return(C)


#inv(L) inverts braids in our setting. The sequence of relatives integers L=[l_1, l_2, ..., l_k] (coresponding to the braid sigma_{|l_1|}^{sign(l_1)}sigma_{|l_2|}^{sign(l_2)}...sigma_{|l_k|)^{sign(l_k)} is sent to the sequence of relatives integers [-l_k, -l_{k-1}, ... -l_1] (coresponding to the braid sigma_{|l_k|)^{sign(-l_k)}sigma_{|l_{k-1}|)^{-sign(l_{k-1})}...sigma_{|l_1|}^{-sign(l_1)}.
def inv(L):
    l=len(L)
    return([-L[l-i-1] for i in range(0,l)])

# commu(U,V) takes two sequences (corrseponding to braids) and return a sequence (corresponding to the commutator of the two braid).
def commu(U,V):
    return(U+V+inv(U)+inv(V))

#function to simplify the sigma_i.sigma_i^{-1]
def simplify(B):
    i=0
    while i<len(B)-1:
        if B[i]==-B[i+1]:
            B.pop(i+1)
            B.pop(i)
            i-=2
        i+=1

#we define the claspers (i,j)
def claspgene(i,j):
    return([j-k for k in range(1,j-i)]+[i,i]+[-i-k for k in range(1,j-i)])

# clasp(I) creates the braid comutator corresponding to the clasper (i_1,...,i_n)=[[A_{i_1,i_n},A_{i_2,i_n}],...]A_{i_{n-1},i_n}]
def clasp(I):
    l=len(I)-1
    T=claspgene(I[0],I[l])
    for i in I[1:l]:
        T=commu(T,claspgene(i,I[l]))
    return(T)

#test verifies if the clasper U of the form (i_1,...,i_n) is in W
def test(U,W):
    for V in W:
        if U==V:
            return(1) 
    return(0)

def conj(T):
    S=[0]+T[:-1]
    return([j+T[-1]-T[-2] for j in S])

    return 
def filtre(C):
    i=0
    while i<len(C):
        T=sorted(C[i][1:])
        W=[]
        for j in T[:-1]:
            T=conj(T)
            W=W+[T]
        m=1
        j=i+1
        while j<len(C):
            if test(sorted(C[j][1:]),W)==1:  
                C.pop(j)
                j-=1
                m+=1
            j+=1
        i+=1

    

def torsion_opti(n):
    L=['lambda']
    B=[-i for i in range(1,n)]
    for p in range(2,n):
        U=B*n
        simplify(U)
        #print('La tresse est de longueur:',len(U))
        C=proj_calcrpz(p,U,[[1,n]])[1:]
        #print('longueur de C avant filtre:',len(C),C)
        filtre(C)
        #print('longueur de C après filtre:',len(C),C)
        for T in C:
            L=[T]+L 
            if T[0]>0:
                B=T[0]*clasp(T[1:])+B
                #print('pour le',i,'ieme générateur on a',len(C),'éléments')
            if T[0]<0:
                B=-T[0]*inv(clasp(T[1:]))+B   
        #print('Les claspeurs de degré',p-1,'à compenser sont donné par:',C)
    U=B*n
    simplify(U)
    #print('La tresse est de longueur:',len(U))
    C=proj_calcrpz(n,U,[[1,n]])[1:]
    return(L,C)